Real-space nanoimaging of THz polaritons in the topological insulator Bi2Se3

Plasmon polaritons in topological insulators attract attention from a fundamental perspective and for potential THz photonic applications. Although polaritons have been observed by THz far-field spectroscopy on topological insulator microstructures, real-space imaging of propagating THz polaritons has been elusive so far. Here, we show spectroscopic THz near-field images of thin Bi2Se3 layers (prototypical topological insulators) revealing polaritons with up to 12 times increased momenta as compared to photons of the same energy and decay times of about 0.48 ps, yet short propagation lengths. From the images we determine and analyze the polariton dispersion, showing that the polaritons can be explained by the coupling of THz radiation to various combinations of Dirac and massive carriers at the Bi2Se3 surfaces, massive bulk carriers and optical phonons. Our work provides critical insights into the nature of THz polaritons in topological insulators and establishes instrumentation and methodology for imaging of THz polaritons.

P lasmon polaritons in metals, doped semiconductors and two-dimensional (2D) materials have wide application potential for field-enhanced spectroscopies, sensing, imaging, and photodetection [1][2][3][4][5][6][7] . Recently, topological insulators (TIs) have been attracting large attention as an alternative class of plasmonic materials, as they can support plasmon polaritons that are formed not only by massive but also by Dirac carriers 8,9 . Dirac plasmon polaritons (DPPs) are electromagnetic modes that can be formed when the massless Dirac carriers at the surfaces of a TI collectively couple to electromagnetic radiation 8,[10][11][12][13] . Due to the 2D nature of these collective excitations, the polariton momentumand thus the field confinementis much larger than that of free space photons of the same energy 8,9 , similar to plasmons in 2D materials such as graphene 14,15 . In addition, spin-momentum locking of the electrons in the TI surface states promises additional unique phenomena, such as spin-polarized plasmon waves 16 . For these reasons, DPPs in TIs have attracted significant interest from both a fundamental and applied perspective [17][18][19][20][21] .
DPPs have been reported experimentally by terahertz (THz) far-field spectroscopy of TI microresonator structures 8,9,22 . Due to the unavoidable presence of massive bulk carriers in TI thin films and crystals 9,13,23,24 , however, the analysis and interpretation of the observed THz resonances has been challenging and controversial. The presence of THz bulk phonon polaritons in the prototypical TI Bi 2 Se 3 further complicates the observation of DPPs 24,25 . Although far-field spectroscopy of TI resonators has provided various fundamental insights, it does not allow for imaging of polariton propagation or mode profiles. Imaging of polaritonsoften performed by scattering-type scanning nearfield optical microscopy (s-SNOM) 26,27 -has proven to be of great importance in the infrared spectral range to distinguish between propagating and localized modes in thin layers and resonator structures, for measuring polariton propagation lengths, phase and group velocities, lifetimes and modal field distributions 6,7,14,15,[27][28][29][30] . However, due to the lack of THz nearfield imaging instrumentation offering high spatial and spectral resolution, as well as a large signal-to-noise ratio, the real-space imaging of THz polaritons is still a challenging task [31][32][33] .
Here, we demonstrate that THz polaritons in TIs can be imaged spectroscopically by s-SNOM employing the tunable monochromatic radiation from a powerful THz gas laser and interferometric detection. Specifically, we performed THz polariton interferometry on epitaxially grown Bi 2 Se 3 films of different thicknesses d. Challenged by the short polariton propagation lengths, we determine the polariton wavevector (and thus dispersion) by complex-valued near-field analysis of our experimental data. Further, using an analytical model, we show that the experimental polariton dispersion can be reproduced when Dirac and massive carriers at the surfaces, massive bulk carriers and optical phonon are taken into account. From propagation length measurements and group velocities determined from the experimental polariton dispersion, we finally determine the decay times of the THz polaritons, amounting to~0.48 ps and thus being comparable or even better than that of typical plasmon decay times in standard (non-encapsulated) graphene.

Results
THz nanoimaging. For real space imaging of polaritons in thin Bi 2 Se 3 films grown by molecular beam epitaxy (Methods section) on sapphire (Al 2 O 3 ), we used a THz s-SNOM (based on a commercial setup from Neaspec, Germany; sketched in Fig. 1a; for details see Supplementary Note 1 and Supplementary Fig. 1), where a metallized atomic force microscope (AFM) tip acts as a THz near-field probe. The tip is illuminated with monochromatic THz radiation from a gas laser (SIFIR-50, Coherent Inc., USA), which is focused with a parabolic mirror. Via the lightning rod effect, the tip concentrates the THz radiation into a nanoscale near-field spot at the tip apex 34 . The momenta of the near fields are large enough to launch polaritons in Bi 2 Se 3 . The tip-launched polaritons that propagate to the edge are reflected at the edge and propagate back to the tip (illustrated in Fig. 1b by the red sine waves). Consequently, by recording the tip-scattered field as function of tip position, we map the interference of forward-and backward-propagating polaritons. Collection and detection of the tip-scattered field is done with the same parabolic mirror and a GaAs-based Schottky diode (WR-0.4ZBD, Virgina Diodes Inc. USA). To obtain background-free near-field signals, the tip is oscillated at a frequency Ω (tapping mode) and the detector signal is demodulated at higher harmonics n of the oscillation frequency, nΩ. Demodulated near-field amplitude and phase signals, s n and φ n , were obtained by synthetic optical holography (SOH), which is based on a Michelson interferometer where the reference mirror (mounted on a delay stage) is translated at a constant velocity along the reference beam path 35,36 . The interferometric detection is key to improve the background suppression and to enable a complex-valued analysis of near-field profiles, which is critical for a reliable measurement of the wavelength of polaritons with short propagation lengths. To increase the signal-to-noise ratio, we used commercial gold tips with a large apex radius 37 of 500 nm (Team Nanotec LRCH). They were oscillated at a frequency of about Ω ≈ 300 kHz with an amplitude of~200 nm.
Representative THz near-field amplitude and phase images, s 3 and φ 3 , of a d = 25-nm-thick Bi 2 Se 3 film are shown in Fig. 1d. The phase image reveals a dark fringe on the Bi 2 Se 3 , which is oriented parallel to the film edge (obtained by scratching the film) and resembles s-SNOM images of short-range plasmon and phonon polaritons observed at mid-IR frequencies on graphene and h-BN, respectively 14,38,39 . In contrast, the simultaneously recorded topography image (Fig. 1c) reveals a homogenous thickness of the Bi 2 Se 3 film, from which we can exclude that the fringe in the THz image is caused by a thickness-dependent dielectric material contrast.
To verify that the dark fringe in Fig. 1d can be attributed to polaritons, we recorded near-field amplitude and phase line profiles perpendicular to the Bi 2 Se 3 edge, s 3 ðxÞ and φ 3 ðxÞ, at different THz frequencies ω. The phase profiles φ 3 ðxÞ are shown in Fig. 1e. We find that the minimum of the near-field phase signal (marked by an arrow, corresponding to the dark fringe of Fig. 1d) shifts towards the Bi 2 Se 3 edge with increasing frequency, supporting our assumption that the near-field signal reveals polaritons of several micrometer wavelength. However, the lack of multiple signal oscillators prevents a straightforward measurement of the polariton wavelength.
Complex-valued analysis of THz line profiles. To establish a procedure for measuring the polariton wavelengths λ p and corresponding wavevector k 0 p ¼ 2π=λ p , we performed a complexvalued analysis of the THz near-field amplitude and phase line profiles, as illustrated in Fig. 2 with data obtained on a 25-nmthick Bi 2 Se 3 film that were recorded at 2.52 THz. We first constructed complex-valued line profiles σ 3 x ð Þ ¼ s 3 x ð Þe iφ 3 x ð Þ and plotted the corresponding trajectories in the complex plane, i.e. as a polar plot where the polar amplitude and phase represent the near-field amplitude s 3 x ð Þ and the near-field phase φ 3 ðxÞ, respectively. We find that the near-field signal describes a spiral (red data in Fig. 2e, based on the line profiles shown in Fig. 2c) around a complex-valued offset C that corresponds to the nearfield signal at large tip-edge distances x. The spiral stems from a harmonic oscillation (describing a circle) whose amplitude decays f Simulation of s-SNOM line profiles: A vertically orientated dipole source (mimicking the tip) is located 1.5 μm above a sheet of conductivity σ (blue, mimicking the Bi 2 Se 3 layer) on Al 2 O 3 . The electric field E z below the dipole at height z NF = 200 nm above the conductivity sheet is calculated and plotted as function of the distance x between the dipole and the sheet edge. g-i Simulated amplitude and phase line profiles analogous to panels c-e. The conductivity was obtained from Eq. (2) with k p = 0.55 + 0.17i μm −1 . For better comparison between the experimental and simulation results, the offset C in the simulations (red data) was replaced by the experimental offset and the phases in panels d and h were set to zero a x = μm. with increasing x, indicating a single propagating mode that is strongly damped. Indeed, after removing the offset (blue data in Fig. 2e), we obtain a monotonically decaying amplitude and a linearly increasing phase signal (Fig. 2d). To verify that the spiral reveals a damped propagating wave, we fitted the complex-valued experimental line profile (red data in Fig. 2e) by which describes the electric field of a back-reflected, radially (i.e. tip-launched) propagating damped wave (black curve in Fig. 2e).
The fitting parameters are A, k p and C. k p is the complex-valued polariton wavevector k p ¼ k 0 p þ ik 00 p ; where k 0 p ¼ 2π=λ p and 1/k 00 p is the propagation length. The complex-valued offset C corresponds to the tip-sample near-field interaction in absence of polaritons that are back-reflected from the Bi 2 Se 3 edge, i.e. when the tip is far away from the edge. A is a complex-valued factor. This offset is present in all polariton maps obtained by s-SNOM and described for example in refs. 40,41 . It can be also seen in our numerical simulations discussed in Fig. 2f-i. After removing the offset C, we obtain Ae À2k 00 , where the term Ae À2k 00 p x = ffiffiffiffiffi 2x p describes a decaying amplitude and the term e i4πx=λ p a linearly increasing phase φ ¼ 4πx=λ p when the distance x to the Bi 2 Se 3 edge increases. The linear relation between distance x and phase φ thus reveals directly the polariton wavelength according to λ p ¼ 4πΔx=Δφ. The fits (black curves in Fig. 2c-e) match well the experimental data, in particular the linear increase of the phase (Fig. 2d) when C is removed, which is the key characteristics of a propagating mode. The fit yields a wavevector of k p = 0.55 + 0.17i μm −1 , corresponding to a normalized wavevector q = k p =k 0 = 10.4 + 3.2i, where k 0 is the photon wavevector. Note that for fitting we excluded the first 200 nm from the edge, in order to avoid a potential influence of tip-edge near-field interaction and edge modes (see discussion below). In the Supplementary Figs. 2 and 3 of the Supplementary Note 2 we show all recorded line profiles and fittings reported in this work.
To verify our analysis of the experimental s-SNOM profiles and the determination of the polariton wavevector k p , we performed well-established numerical model simulations 30,42 . As illustrated in Fig. 2f, the s-SNOM tip is described by a vertically orientated dipole source and the Bi 2 Se 3 layer by a 2D sheet of an optical conductivity σ (blue layer in Fig. 2f). The electric field E z ¼ E z e iφ z (describing the s-SNOM signal) below the dipole is calculated and plotted as function of the distance x between the dipole source and the sheet edge (mimicking the scanning of the tip, Fig. 2g). To obtain the sheet conductivity σ, we assume that optical polariton modes are probed in our experiment (where the polariton fields normal to the film have opposite sign at the top and bottom surface; for further discussion see below). In this case, σ can be obtained from the dispersion relation of polaritons in a 2D sheet within the large momentum approximation 43,44 : where q is the normalized complex-valued polariton wavevector along the film, ω the frequency, and ε sub = 10 the permittivity of the Al 2 O 3 substrate 9,24,45 (see Supplementary Note 3 and Supplementary Fig. 4). We note that the approximation of the layer by a 2D sheet conductivity is justified, despite Bi 2 Se 3 being an anisotropic material hosting hyperbolic polaritons, as the layers are much thinner than the polariton wavelength and the damping of the polaritons is rather high (for further details see Supplementary Fig. 5 and Supplementary Notes 3 and 4). For k p = 0.55 + 0.17i μm −1 (according to our analysis of the experimental s-SNOM line profiles in Fig. 2c-e), we obtain the simulated near-field line profiles shown in Fig. 2g-i. An excellent agreement with the experimental s-SNOM line profiles is found. Particularly, complex-valued fitting of the simulated line profiles by a radially decaying wave (according to Eq. (1)) yields a polariton wavevector k p = 0.63 + 0.15i μm −1 , which closely matches the value determined from the experimental s-SNOM line profiles. The simulations thus confirm that the experimental s-SNOM line profiles reveal a polariton mode that (i) is launched by the tip, (ii) propagates as a damped wave radially along the Bi 2 Se 3 film, (iii) reflects at the edge of the film back to the tip, and (iv) is scattered by the tip. For a demonstration of our conclusions, we show the electric near-field distribution around the dipole source placed above the conductivity sheet, Re E z x; y À Á Â Ã : When the dipole is placed inside the sheet, i.e. far away from any edge, we clearly observe radially propagating wavefronts (Fig. 3a, upper panel). Most important, fitting of the field distribution along the horizontal dashed red line by Re Ae ik p x = ffiffi ffi x p Â Ã yields k p = 0.60 + 0.15i μm −1 (Fig. 3a,  lower panel), which agrees well with the wavevectors k p obtained from the experimental and simulated s-SNOM line profiles. The slight discrepancies between the various wavevectors (<15%) may be attributed to the excitation of an edge mode when the tip comes into close proximity of the sheet edge. The edge mode can be actually recognized in the simulations when the dipole is located at the sheet edge (Fig. 3b, upper panel). Its wavelength is slightly reduced compared to that of the sheet mode, and its fields are strongly confined to the edge, similarly to what has been observed for plasmon and phonon polariton modes in graphene and h-BN flakes 30,46,47 . Most important, since the edge mode propagates exclusively along the edge and its field is strongly confined to the edge, its contribution in the experimental and simulated s-SNOM line profiles perpendicular to the edge (shown in Fig. 2) is minor, when we allow small uncertainties below 15%. Analysis of the THz polariton dispersion. In Fig. 4a we show the phase images of two Bi 2 Se 3 films of different thicknesses recorded at different frequencies, which were used to determine the polariton dispersions according to the procedure described in Fig. 2 (red symbols in Fig. 4b; q = k p =k 0 is the polariton wavevector k p normalized to the photon wavevector k 0 ). As is typical for polaritons, we find that Re½q increases with increasing frequency ω and with decreasing film thickness d. To understand the physical origin of the polaritons, we compare the experimental results with analytical calculations of the polariton dispersion employing Eq. (2) and various conductivity models describing the Bi 2 Se 3 film (see Supplementary Note 4, Supplementary Fig. 6 and Supplementary Table 1). For modeling of the conductivity, we performed Hall measurements of Bi 2 Se 3 films of different thicknesses (see Supplementary Fig. 7 and Supplementary Notes 5 and 6), yielding an effective 2D carrier concentration of about n 2D;Hall = 2.5 × 10 13 cm −2 for layers with a thickness of~25 nm ( Supplementary Fig. 8).
First, we assume that only optical phonons (OP) and Dirac carriers (DC) located on both film surfaces contribute to the conductivity, yielding σ ¼ σ model OPþDC ¼ ωd 4πi ε phonon þ 2σ Dirac ; where ε phonon is the bulk dielectric function of Bi 2 Se 3 including optical phonons. σ Dirac is the sheet conductivity of one Bi 2 Se 3 surface (see Supplementary Eq. (8)) 11,12,25 , assuming that the sheet carrier concentration at one surface is n Dirac ¼ n 2D;Hall =2 = 1.25 × 10 13 cm −2 (independent of the thickness). The resulting dispersions are shown by the orange curves   in Fig. 4b. We find that the calculated wavevectors are significantly larger than the experimental values, from which we conclude that pure Dirac plasmon polaritons coupled to phonon polaritons cannot explain the experimental dispersion. As a second case, we assume that all carriers are bulk carriers (BC), yielding σ ¼ σ model OPþBC ¼ ωd 4πi ðε phonon þ ε Drude Þ; where ε Drude is the Drude contribution to the bulk dielectric function (see Supplementary Note 4). In this case, we use an effective three-dimensional (3D) concentration of the massive carriers according to n bulk = n 2D;Hall /25 nm = 1 × 10 19 cm −3 . We obtain the dispersions shown by the light blue curves in Fig. 4b. For the 25 nm thick film a reasonable match of the experimental polariton wavevectors is found, however, not for the 60 nm thick film, revealing that polaritons comprising only bulk carriers (BC) and optical phonons (OP) cannot explain the polariton dispersions either. A similar observation is made (green curves in Fig. 4b) when we assume that all carriers stem from a massive 2D electron gas (2DEG) -which is known to exist in TIs due to surface band bending 18,48,49 11,12,49 is the sheet conductivity of one Bi 2 Se 3 surface with n 2DEG ¼ n 2D;Hall =2 = 1.25 × 10 13 cm −2 .
We next assume that both massive bulk carriers and surface carriers (either Dirac or massive 2DEG carriers) contribute to the conductivity, σ ¼ σ model , respectively. Note that massive bulk carriers in thin films are barely captured by the Hall measurements (due to their supposedly smaller mobility compared to that of the Dirac carriers; see Supplementary Note 6). We thus assign the total Hall-measured 2D concentration (n 2D;Hall = 2.5 × 10 13 cm −2 ) fully to either Dirac or massive 2DEG carriers for both the 25 nm and 60 nm thick film, and consider an additional massive bulk carrier concentration n bulk . From Hall measurements of thick Bi 2 Se 3 films we estimate n bulk = 2.15 × 10 18 cm −3 (Supplementary Note 6), yielding the bulk Drude contribution ε Drude . The calculated dispersions (black and purple lines Fig. 4c, labeled OP+BC+2DEG and OP+BC+DC, respectively) again do not match the experimental results (red symbols).
In the following, we attempt to fit the experimental dispersions by various parameter variations. We first added an additional 2DEG contribution, such that σ ¼ σ model (red lines Fig.  4d). Using for each surface the carrier concentrations n bulk = 2.15 × 10 18 cm −3 and n Dirac ¼ 1.25 × 10 13 cm −2 (from the Hall measurements), we obtain the fitting parameter n 2DEG ¼ 0.375 × 10 13 cm −2 for each surface. We note that in Hall measurements we cannot separate Dirac and massive carriers directly (Supplementary Note 6) to confirm these carrier concentrations, but they are close to the numbers reported in literature 9,23,50-52 . Interestingly, the experimental dispersions can be also fitted without considering a 2DEG (employing the conductivity model σ ¼ σ model OPþBCþDC , light purple line in Fig. 4d). However, we have to assume an increased bulk carrier concentration of n bulk ¼ 3.72 × 10 18 cm −3 (fitting parameter; n Dirac as before). Although the required bulk carrier concentration is nearly twice as high as the one estimated from our Hall measurements, it represents a reasonable value reported in literature 9,23,24,49,50 , which may not be fully revealed by Hall measurements (see discussion in Supplementary Note 6). We also fitted the experimental dispersions without considering Dirac carriers employing the conductivity model σ ¼ σ model OPþBCþ2DEG (gray line in Fig. 4d) with n bulk = 2.15 × 10 18 cm −3 and n 2DEG being the fit parameter. A good matching of the experimental dispersion is achieved for n 2DEG ¼ 0.95 × 10 13 cm −2 for each surface. However, this value of n total;2DEG ¼ 1.9 × 10 13 cm −2 is significantly higher than what has been reported in literature 18,49,51 . Altogether, we conclude from our systematic dispersion analysis that an unambiguous clarification of the nature and concentration of carriers forming the polaritons is difficult without additional experiments where the concentrations of the different carriers can be measured separately. However, such measurements are challenging to carry out at room temperature due to thermal smearing, and low-temperature measurements are unlikely to be accurate at room temperature due to thermal excitations. On the other hand, considering that Bi 2 Se 3 growth is highly reproducible and that Dirac carriers have been verified in samples like ours 53 , these carriers may contribute to the signal.
Polariton propagation length and lifetime. From the complexvalued fitting of the s-SNOM line profiles we also obtain the propagation length of the polaritons, L ¼ 1=k 00 p . For the 25 nm thick film we find L = 6 μm and accordingly the amplitude decay time τ ¼ L=v g = 0.48 ps, which is similar to the decay times measured by far-field extinction spectroscopy of Bi 2 Se 3 ribbons 9 . The group velocities v g were obtained from the polariton dispersion according to v g ¼ dω=dk 0 p = 0.042c. Interestingly, the polariton decay time in Bi 2 Se 3 is comparable or even larger than that of graphene plasmons at infrared frequencies 14,15,28 . On the other hand, the inverse damping ratio of the Bi 2 Se 3 polaritons, γ À1 = k 0 p =k 00 p = 3.2, and their relative propagation length, L=λ p ¼ 1 2πγ ¼ 0.5, is significantly smaller than that of the infrared graphene plasmons (γ À1 = 5) 15 . To understand the small relative propagation length of the Bi 2 Se 3 polaritons, we express it as a function of the polariton wavevector k p ¼ qω=c, group velocity v g and decay time τ: From Eq. (3) it becomes clear that for a given τ, v g and q, the relative propagation length decreases with decreasing frequency; that is simply because the temporal oscillation period becomes longer. Since typical THz frequencies are more than one order of magnitude smaller than infrared frequencies, one can expect, generally, that the relative propagation length of THz polaritons in thin layers (including 2D materials) is significantly smaller than that of infrared polaritons. For an illustration, we show in Fig. 5a the calculated relative propagation length L=λ p of 2.52 THz polaritons of 0.48 ps amplitude decay time as function of Re½q and v g . The white symbol marks the relative propagation length of the THz polaritons observed in the 25-nm-thick Bi 2 Se 3 film. We find that propagation lengths of more than wavelength, L=λ p > 1, are possible only for large group velocities (>0.1c) when the normalized polariton wavevector (i.e. polariton confinement) is moderate (Re½q > 10). To achieve L=λ p > 1 for group velocities below 0.05c, large polariton wavevectors with Re½q > 20 are required. For comparison, we show in Fig. 5b the relative propagation length for polaritons of 0.6 ps amplitude decay time. Plasmon polaritons of such rather exceptionally large amplitude decay time were observed experimentally in high-quality graphene encapsulated in h-BN (marked by black symbol). Only because of their very large confinement (Re½q = 70; owing to coupling with adjacent metallic gate electrodes), these polaritons possess a large relative propagation length of L=λ p ¼ 1.5, although their group velocity is rather small (v g = 0.014c). Generally, we conclude from Eq. (3) and Fig. 5 that the relative propagation lengths of THz polaritons are generally short, unless THz polaritons of extraordinary long decay times 33 , large wavevectors or large group velocities are studied.

Discussion
We note that in our experiments we could only observe the optical polariton modes, although s-SNOM in principle can map acoustic polariton modes as well 31 . We explain the absence of acoustic polariton modes (where the sign of the effective surface charges is opposite on both surfaces) by their extremely short wavelengths, which may prevent efficient coupling with the probing tip. Further, the acoustic modes might be strongly overdamped due to the presence of bulk carriers. In the future, sharper tips and TIs of lower bulk carrier concentration which can be grown using a buffer layer technique 54 and may allow the study of the ultra-confined acoustic modes in real space.
In summary, we demonstrated an instrumentation for s-SNOM that allows for spectroscopic nanoimaging of thin-film polaritons around 2.5 THz, even in case of weak polaritonic image contrasts. We applied it to record real-space images of THz polaritons in the TI Bi 2 Se 3 . Despite the short polariton propagation, we could measure the polariton dispersion and propagation length, owing to complex-valued analysis of near-field line profiles. In the future, the highly specific signature of polaritonic spatial signal oscillations in the complex planerepresenting a spiralcould be also applied to distinguish them from nonpolaritonic spatial signal oscillations that, for example, are caused by spatial variations of dielectric material properties or by laser intensity fluctuations. Using dispersion calculations based on various optical conductivity models, we found that the polaritons can be explained by simultaneous coupling of THz radiation to various combination of Dirac carriers, massive 2DEG carriers, massive bulk carriers and optical phonons. During the revision of our manuscript, another THz s-SNOM study of polaritons in TIs was published 55 , reporting that massive 2DEGs need to be considered when interpreting THz polaritons in Bi 2 Se 3 . We note, however, that the contribution of massive carriers may be strongly reduced or even absent 2,13 by growing the Bi 2 Se 3 by alternative methods, for example by using a buffer layer technique 52,56 . Beyond s-SNOM-based dispersion analysis as demonstrated here, our work paves the way for studying THz polaritons on other TI materials, 2D materials or 2DEGs, such as the mapping of modal field patterns in resonator structures 30 and moiré superlattices 57 , or the directional propagation on in-plane anisotropic natural materials 42 and metasurfaces 58 .

Methods
Sample preparation. Films of Bi 2 Se 3 are grown via molecular beam epitaxy (Veeco GenXplor MBE system) on single-side polished sapphire substrate (0001) plane (10 mm × 10 mm × 0.5 mm, MTI Corp., U.S.A.). All films are grown at the same substrate temperature as measured by a non-contact thermocouple (325°C), growth rate (0.8 nm/min), and selenium: bismuth flux ratio as measured by an ion gauge (≈50), and the selenium cell has a high-temperature cracker zone set at 900°C 52 . Film thicknesses are determined via x-ray reflection (XRR) measurement. X-ray diffraction further confirmed that only a single phase and one orientation (0001) of Bi 2 Se 3 has been epitaxially grown on c-direction on the substrate. The sheet concentration (~3.0 × 10 13 cm −2 ) for 120-nm-thick Bi 2 Se 3 are obtained via Hall effect measurement in a van der Pauw configuration at room temperature (see Supplementary Note 5), with error bars~6%.

Data availability
Data that support the results of this work are available upon reasonable request from the corresponding author.